##########################################################################
# June 20 2024
# Replication File
# Rebecca Cordell and Kashmiri Medhi
# Transnational Repression: International Cooperation in Silencing Dissent
# International Studies Quarterly
# https://academic.oup.com/isq
##########################################################################

# Clear work environment
rm(list=ls())

# Install Packages
#install.packages("stats")
#install.packages("countrycode")
#install.packages("RColorBrewer")
#install.packages("rworldmap")
#install.packages("pscl")
#install.packages("texreg")
#install.packages("lmtest")
#install.packages("ggeffects")
#install.packages("ggplot2")
#install.packages("gamlss")

# Required Packages
require(stats)
require(countrycode)
require(RColorBrewer)
require(rworldmap)
require(pscl)
require(texreg)
require(lmtest)
require(ggeffects)
require(ggplot2)
require(gamlss)

options(scipen=999)

###############
# Read in data
###############
tr<-read.csv("CordellMedhi_2024_TransnationalRepression.csv", header= TRUE, stringsAsFactors = FALSE)

# Reduce to sample of host countries admitting refugees or asylum seekers from the origin country that year
tr_rd<-tr[tr$refugees!=0,]

###################################################################
# Figure 1: Origin Countries, Frequency of Transnational Repression
###################################################################

# Reduce to origin countries and host countries who have engaged in transnational repression according to Freedom House (2021) Transnational Repression Database
tr_all<-tr[tr$tr_count_all!=0,]

# Reduce to origin countries and count the total number of instances
fig1_df<-tr_all[, c("iso3c2","tr_count_all")]
fig1_df<-aggregate(fig1_df$tr_count_all, by=list(Category=fig1_df$iso3c2), FUN=sum)

# Join data and country map
map.data <- joinCountryData2Map(        
  fig1_df,                         
  joinCode = "ISO3",                    
  nameJoinColumn = "Category",           
  verbose = TRUE)
map.data <- subset(map.data, continent != "Antarctica")

# Create colour palette
palette = colorRampPalette(brewer.pal(n=7, name='YlOrRd'))(7)

# Print Figure 1
mapCountryData(map.data,                                     
               nameColumnToPlot="x",                 
               addLegend=TRUE,  
               colourPalette=palette,
               catMethod="fixedWidth",                                                 
               missingCountryCol="white", 
               mapTitle="", 
               borderCol="black")

#################################################################
# Figure 2: Host Countries, Frequency of Transnational Repression
#################################################################

# Reduce to origin countries and count the total number of instances
fig2_df<-tr_all[, c("iso3c1","tr_count_all")]
fig2_df<-aggregate(fig2_df$tr_count_all, by=list(Category=fig2_df$iso3c1), FUN=sum)

# Join data and country map
map.data <- joinCountryData2Map(        
  fig2_df,                         
  joinCode = "ISO3",                    
  nameJoinColumn = "Category",           
  verbose = TRUE)
map.data <- subset(map.data, continent != "Antarctica")

# Create colour palette
palette = colorRampPalette(brewer.pal(n=7, name='YlOrRd'))(7)

# Print Figure 2
mapCountryData(map.data,                                     
               nameColumnToPlot="x",                 
               addLegend=TRUE,  
               colourPalette=palette,
               catMethod="fixedWidth",                                                 
               missingCountryCol="white", 
               mapTitle="", 
               borderCol="black")

###################################################################################################
# Table 1: Zero-inflated negative binomial regression, count of transnational repression incidences
###################################################################################################

# Model 1: Baseline model, independent variables
tab1_m1 <- zeroinfl(tr_count_det_ren_dep ~ rol_inverse + trade_total_log | 1, dist = 'negbin', data = tr_rd)

# Model 2: Full model, independent variables and control variables
tab1_m2 <- zeroinfl(tr_count_det_ren_dep ~ rol_inverse + trade_total_log + alliance + refugees_log + min_dist_log + democracy + state_capac + gdppc_log + pop_log + right_orient + civ_lib + orign_country_past + press_freedom | 1, dist = 'negbin', data = tr_rd)

# Print Regression Table
screenreg(list(tab1_m1, tab1_m2), custom.model.names = c("Model 1", "Model 2"), stars = c(0.01, 0.05, 0.1), digits = 3)

# Observations
tab1_n_m1<-tab1_m1$n
tab1_n_m2<-tab1_m2$n

# Log Likelihood
tab1_loglik_m1<-tab1_m1$loglik
tab1_loglik_m2<-tab1_m2$loglik

# theta 
tab1_theta_m1<-tab1_m1$theta
tab1_theta_m2<-tab1_m2$theta

# Prob. > X2
tab1_px2_m1<-lrtest(tab1_m1)
tab1_px2_m1<-tab1_px2_m1$`Pr(>Chisq)`[2]
tab1_px2_m2<-lrtest(tab1_m2)
tab1_px2_m2<-tab1_px2_m2$`Pr(>Chisq)`[2]

# Create Model Fit Stats Table
tab1_m1_fit<-rbind(tab1_n_m1, tab1_loglik_m1, tab1_theta_m1, tab1_px2_m1)
tab1_m2_fit<-rbind(tab1_n_m2, tab1_loglik_m2, tab1_theta_m2, tab1_px2_m2)
tab1_fit<-round(cbind(tab1_m1_fit,tab1_m2_fit), digits=3)

# Print Model Fit Stats
tab1_fit

##############################################################################
# Figure 3: Rule of Law (inverse), Predicted Count of Transnational Repression
##############################################################################

# Predicted Count, Rule of Law (inverse), Model 2
rol_pred <- ggpredict(tab1_m2, c("rol_inverse"), full.data = TRUE)

# "In Figure 3, moving from the minimum to the maximum value of Rule of Law (inverse) variable increases the predicted count of transnational repression from 0 (countries with a strong rule of law) to 1.52 (countries with a weak rule of law)." 
round(rol_pred$predicted, digits=2)

# Print Figure 3
fig3<-ggplot(data=rol_pred, aes(x=x, y=predicted, group=1)) +
  geom_line()+
  geom_point() + theme_bw() + labs(y = "Transnational Repression", x = "Rule of Law (inverse)") +
  geom_ribbon(aes(ymin = conf.low, ymax = conf.high), alpha = .1)
fig3

##########################################################################
# Figure 4: Total Trade (log), Predicted Count of Transnational Repression
##########################################################################

# Predicted Count, Total Trade (log), Model 2
trade_pred <- ggpredict(tab1_m2, c("trade_total_log"), full.data = TRUE)

# "In Figure 4, shifting the volume of dyadic trade from the minimum value to the maximum value in our data increases the predicted count of transnational repression from 0 (low levels of trade) to 1.22 (high levels of trade)."
round(trade_pred$predicted, digits=2)

# Print Figure 4
fig4<-ggplot(data=trade_pred, aes(x=x, y=predicted)) +
  geom_line()+
  geom_point() + theme_bw() + labs(y = "Transnational Repression", x = "Total Trade (log)") +
  geom_ribbon(aes(ymin = conf.low, ymax = conf.high), alpha = .1)
fig4

############################################################################################################################
# Figure 5: Interaction between the Rule of Law (inverse) and Total Trade (log), Predicted Count of Transnational Repression
############################################################################################################################

# Model 1: Baseline model, independent variables with interaction
fig5_m1 <- zeroinfl(tr_count_det_ren_dep ~ rol_inverse*trade_total_log + rol_inverse + trade_total_log | 1, dist = 'negbin', data = tr_rd)

# Model 2: Full model, independent variables with interaction and control variables
fig5_m2 <- zeroinfl(tr_count_det_ren_dep ~ rol_inverse*trade_total_log + rol_inverse + trade_total_log + alliance + refugees_log + min_dist_log + democracy + state_capac + gdppc_log + pop_log + right_orient + civ_lib + orign_country_past + press_freedom | 1, dist = 'negbin', data = tr_rd)

# Predicted Count, Rule of Law (inverse)*Total Trade (log), Model 2
interact_pred<-ggpredict(fig5_m2, c("rol_inverse","trade_total_log [meansd]"))
interact_pred$group<-as.character(interact_pred$group)
interact_pred$group[interact_pred$group=="13.42"]<-"-1 sd"
interact_pred<-interact_pred[interact_pred$group!="18.01",] #mean
interact_pred$group[interact_pred$group=="22.6"]<-"+1 sd"
interact_pred$group <- factor(interact_pred$group, levels = c("-1 sd", "+1 sd"))

# "However, when host countries have high economic ties to origin countries and a weak rule of law, the predicted count of transnational repression dramatically increases from 0 (no instances) to 34 occurrences."
round(interact_pred$predicted)

# Print Figure 5
fig5 <- ggplot(interact_pred, aes(x = x, y = predicted, group = group)) +
  geom_line(aes(col = group)) +
  geom_point() +
  geom_ribbon(aes(ymin = conf.low, ymax = conf.high), alpha = 0.1) +
  theme_bw() +
  labs(y = "Transnational Repression", x = "Rule of Law (inverse)", color = "Total Trade (log)") +
  guides(colour = guide_legend(reverse=T))
fig5

##################################################################################################################
# Table A.1: Zero-inflated negative binomial regression, count of transnational repression incidences, full sample
##################################################################################################################

# Model 1: Baseline model, independent variables
taba1_m1 <- zeroinfl(tr_count_det_ren_dep ~ rol_inverse + trade_total_log | 1, dist = 'negbin', data = tr)

# Model 2: Full model, independent variables and control variables
taba1_m2 <- zeroinfl(tr_count_det_ren_dep ~ rol_inverse + trade_total_log + alliance + refugees_log + min_dist_log + democracy + state_capac + gdppc_log + pop_log + right_orient + civ_lib + orign_country_past + press_freedom | 1, dist = 'negbin', data = tr)

# Print Regression Table
screenreg(list(taba1_m1, taba1_m2), custom.model.names = c("Model 1", "Model 2"), stars = c(0.01, 0.05, 0.1), digits = 3)

# Observations
taba1_n_m1<-taba1_m1$n
taba1_n_m2<-taba1_m2$n

# Log Likelihood
taba1_loglik_m1<-taba1_m1$loglik
taba1_loglik_m2<-taba1_m2$loglik

# theta 
taba1_theta_m1<-taba1_m1$theta
taba1_theta_m2<-taba1_m2$theta

# Prob. > X2
taba1_px2_m1<-lrtest(taba1_m1)
taba1_px2_m1<-taba1_px2_m1$`Pr(>Chisq)`[2]
taba1_px2_m2<-lrtest(taba1_m2)
taba1_px2_m2<-taba1_px2_m2$`Pr(>Chisq)`[2]

# Create Model Fit Stats Table
taba1_m1_fit<-rbind(taba1_n_m1, taba1_loglik_m1, taba1_theta_m1, taba1_px2_m1)
taba1_m2_fit<-rbind(taba1_n_m2, taba1_loglik_m2, taba1_theta_m2, taba1_px2_m2)
taba1_fit<-round(cbind(taba1_m1_fit,taba1_m2_fit), digits=3)

# Print Model Fit Stats
taba1_fit

#########################################################################################################################################
# Table A.2: Zero-inflated negative binomial regression, count of transnational repression incidences, higher annual refugee flows sample
#########################################################################################################################################

# Reduce to sample of host countries with higher annual refugees flows (above the median value)
summary(tr_rd$refugees)
tr_rd_high<-tr_rd[tr_rd$refugees>=30,]

# Model 1: Baseline model, independent variables
taba2_m1 <- zeroinfl(tr_count_det_ren_dep ~ rol_inverse + trade_total_log | 1, dist = 'negbin', data = tr_rd_high)

# Model 2: Full model, independent variables and control variables
taba2_m2 <- zeroinfl(tr_count_det_ren_dep ~ rol_inverse + trade_total_log + alliance + refugees_log + min_dist_log + democracy + state_capac + gdppc_log + pop_log + right_orient + civ_lib + orign_country_past + press_freedom | 1, dist = 'negbin', data = tr_rd_high)

# Print Regression Table
screenreg(list(taba2_m1, taba2_m2), custom.model.names = c("Model 1", "Model 2"), stars = c(0.01, 0.05, 0.1), digits = 3)

# Observations
taba2_n_m1<-taba2_m1$n
taba2_n_m2<-taba2_m2$n

# Log Likelihood
taba2_loglik_m1<-taba2_m1$loglik
taba2_loglik_m2<-taba2_m2$loglik

# theta 
taba2_theta_m1<-taba2_m1$theta
taba2_theta_m2<-taba2_m2$theta

# Prob. > X2
taba2_px2_m1<-lrtest(taba2_m1)
taba2_px2_m1<-taba2_px2_m1$`Pr(>Chisq)`[2]
taba2_px2_m2<-lrtest(taba2_m2)
taba2_px2_m2<-taba2_px2_m2$`Pr(>Chisq)`[2]

# Create Model Fit Stats Table
taba2_m1_fit<-rbind(taba2_n_m1, taba2_loglik_m1, taba2_theta_m1, taba2_px2_m1)
taba2_m2_fit<-rbind(taba2_n_m2, taba2_loglik_m2, taba2_theta_m2, taba2_px2_m2)
taba2_fit<-round(cbind(taba2_m1_fit,taba2_m2_fit), digits=3)

# Print Model Fit Stats
taba2_fit

########################################################################################
# Table A.3: Negative binomial regression, count of transnational repression incidences
#######################################################################################

# Reduce to sample of complete cases
taba3_m1_df<-tr_rd[c("tr_count_det_ren_dep","trade_total_log","rol_inverse")]
taba3_m1_df <- taba3_m1_df[complete.cases(taba3_m1_df),]

taba3_m2_df<-tr_rd[c("tr_count_det_ren_dep","trade_total_log","rol_inverse","alliance","refugees_log","min_dist_log","democracy","state_capac","gdppc_log","pop_log","right_orient","civ_lib","orign_country_past","press_freedom")]
taba3_m2_df <- taba3_m2_df[complete.cases(taba3_m2_df),]

# Model 1: Baseline model, independent variables
taba3_m1 <- gamlss(tr_count_det_ren_dep ~ rol_inverse + trade_total_log, data=taba3_m1_df, family=NBII)

# Model 2: Full model, independent variables and control variables
taba3_m2 <- gamlss(tr_count_det_ren_dep ~ rol_inverse + trade_total_log + alliance + refugees_log + min_dist_log + democracy + state_capac + gdppc_log + pop_log + right_orient + civ_lib + orign_country_past + press_freedom, data=taba3_m2_df, family=NBII)

# Print Regression Table
screenreg(list(taba3_m1, taba3_m2), custom.model.names = c("Model 1", "Model 2"), stars = c(0.01, 0.05, 0.1), digits = 2)

# Observations
taba3_n_m1<-taba3_m1$n
taba3_n_m2<-taba3_m2$n

# Log Likelihood
taba3_loglik_m1<-lrtest(taba3_m1)
taba3_loglik_m1<-taba3_loglik_m1$`LogLik`[1]
taba3_loglik_m2<-lrtest(taba3_m2)
taba3_loglik_m2<-taba3_loglik_m2$`LogLik`[1]

# Prob. > X2
taba3_px2_m1<-lrtest(taba3_m1)
taba3_px2_m1<-taba3_px2_m1$`Pr(>Chisq)`[2]
taba3_px2_m2<-lrtest(taba3_m2)
taba3_px2_m2<-taba3_px2_m2$`Pr(>Chisq)`[2]

# Create Model Fit Stats Table
taba3_m1_fit<-rbind(taba3_n_m1, taba3_loglik_m1, taba3_px2_m1)
taba3_m2_fit<-rbind(taba3_n_m2, taba3_loglik_m2, taba3_px2_m2)
taba3_fit<-round(cbind(taba3_m1_fit,taba3_m2_fit), digits=3)

# Print Model Fit Stats
taba3_fit

#########################################################################
# Table A.4: Descriptive statistics of independent and control variables
########################################################################

# Remove missing values
trade_total_log <- tr_rd[!is.na(tr_rd$trade_total_log),]
rol_inverse <- tr_rd[!is.na(tr_rd$rol_inverse),]
democracy <- tr_rd[!is.na(tr_rd$democracy),]
gdppc_log <- tr_rd[!is.na(tr_rd$gdppc_log),]
alliance <- tr_rd[!is.na(tr_rd$alliance),]
right_orient <- tr_rd[!is.na(tr_rd$right_orient),]
pop_log <- tr_rd[!is.na(tr_rd$pop_log),]
refugees_log <- tr_rd[!is.na(tr_rd$refugees_log),]
civ_lib <- tr_rd[!is.na(tr_rd$civ_lib),]
orign_country_past <- tr_rd[!is.na(tr_rd$orign_country_past),]
press_freedom <- tr_rd[!is.na(tr_rd$press_freedom),]
state_capac <- tr_rd[!is.na(tr_rd$state_capac),]
min_dist_log <- tr_rd[!is.na(tr_rd$min_dist_log),]

# N
n_trade_total_log<-length(trade_total_log$trade_total_log)
n_rol_inverse<-length(rol_inverse$rol_inverse)
n_democracy<-length(democracy$democracy)
n_gdppc_log<-length(gdppc_log$gdppc_log)
n_alliance<-length(alliance$alliance)
n_right_orient<-length(right_orient$right_orient)
n_pop_log<-length(pop_log$pop_log)
n_refugees_log<-length(refugees_log$refugees_log)
n_civ_lib<-length(civ_lib$civ_lib)
n_orign_country_past<-length(orign_country_past$orign_country_past)
n_press_freedom<-length(press_freedom$press_freedom)
n_state_capac<-length(state_capac$state_capac)
n_min_dist_log<-length(min_dist_log$min_dist_log)

# Mean
mean_trade_total_log<-mean(trade_total_log$trade_total_log) 
mean_rol_inverse<-mean(rol_inverse$rol_inverse)
mean_democracy<-mean(democracy$democracy)
mean_gdppc_log<-mean(gdppc_log$gdppc_log)
mean_alliance<-mean(alliance$alliance)
mean_right_orient<-mean(right_orient$right_orient)
mean_pop_log<-mean(pop_log$pop_log)
mean_refugees_log<-mean(refugees_log$refugees_log)
mean_civ_lib<-mean(civ_lib$civ_lib) 
mean_orign_country_past<-mean(orign_country_past$orign_country_past) 
mean_press_freedom<-mean(press_freedom$press_freedom) 
mean_state_capac<-mean(press_freedom$state_capac) 
mean_min_dist_log<-mean(min_dist_log$min_dist_log) 

# Standard Deviation
sd_trade_total_log<-sd(trade_total_log$trade_total_log) 
sd_rol_inverse<-sd(rol_inverse$rol_inverse)
sd_democracy<-sd(democracy$democracy)
sd_gdppc_log<-sd(gdppc_log$gdppc_log)
sd_alliance<-sd(alliance$alliance)
sd_right_orient<-sd(right_orient$right_orient)
sd_pop_log<-sd(pop_log$pop_log)
sd_refugees_log<-sd(refugees_log$refugees_log)
sd_civ_lib<-sd(civ_lib$civ_lib) 
sd_orign_country_past<-sd(orign_country_past$orign_country_past) 
sd_press_freedom<-sd(press_freedom$press_freedom) 
sd_state_capac<-sd(state_capac$state_capac) 
sd_min_dist_log<-sd(press_freedom$min_dist_log) 

# Min
min_trade_total_log<-min(trade_total_log$trade_total_log) 
min_rol_inverse<-min(rol_inverse$rol_inverse)
min_democracy<-min(democracy$democracy)
min_gdppc_log<-min(gdppc_log$gdppc_log)
min_alliance<-min(alliance$alliance)
min_right_orient<-min(right_orient$right_orient)
min_pop_log<-min(pop_log$pop_log)
min_refugees_log<-min(refugees_log$refugees_log)
min_civ_lib<-min(civ_lib$civ_lib) 
min_orign_country_past<-min(orign_country_past$orign_country_past) 
min_press_freedom<-max(press_freedom$press_freedom) 
min_state_capac<-min(press_freedom$state_capac) 
min_min_dist_log<-min(press_freedom$min_dist_log) 

# Max
max_trade_total_log<-max(trade_total_log$trade_total_log) 
max_rol_inverse<-max(rol_inverse$rol_inverse)
max_democracy<-max(democracy$democracy)
max_gdppc_log<-max(gdppc_log$gdppc_log)
max_alliance<-max(alliance$alliance)
max_right_orient<-max(right_orient$right_orient)
max_pop_log<-max(pop_log$pop_log)
max_refugees_log<-max(refugees_log$refugees_log)
max_civ_lib<-max(civ_lib$civ_lib) 
max_orign_country_past<-max(orign_country_past$orign_country_past) 
max_press_freedom<-max(press_freedom$press_freedom) 
max_state_capac<-min(press_freedom$state_capac) 
max_min_dist_log<-min(press_freedom$min_dist_log) 

# Create table
all_trade_total_log<-cbind(n_trade_total_log, mean_trade_total_log, sd_trade_total_log, min_trade_total_log, max_trade_total_log)
all_rol_inverse<-cbind(n_rol_inverse, mean_rol_inverse, sd_rol_inverse, min_rol_inverse, max_rol_inverse)
all_democracy<-cbind(n_democracy, mean_democracy, sd_democracy, min_democracy, max_democracy)
all_gdppc_log<-cbind(n_gdppc_log, mean_gdppc_log, sd_gdppc_log, min_gdppc_log, max_gdppc_log)
all_alliance<-cbind(n_alliance, mean_alliance, sd_alliance, min_alliance, max_alliance)
all_right_orient<-cbind(n_right_orient, mean_right_orient, sd_right_orient, min_right_orient, max_right_orient)
all_pop_log<-cbind(n_pop_log, mean_pop_log, sd_pop_log, min_pop_log, max_pop_log)
all_refugees_log<-cbind(n_refugees_log, mean_refugees_log, sd_refugees_log, min_refugees_log, max_refugees_log)
all_civ_lib<-cbind(n_civ_lib, mean_civ_lib, sd_civ_lib, min_civ_lib, max_civ_lib)
all_orign_country_past<-cbind(n_orign_country_past, mean_orign_country_past, sd_orign_country_past, min_orign_country_past, max_orign_country_past)
all_press_freedom<-cbind(n_press_freedom, mean_press_freedom, sd_press_freedom, min_press_freedom, max_press_freedom)
all_state_capac<-cbind(n_state_capac, mean_state_capac, sd_state_capac, min_state_capac, max_state_capac)
all_min_dist_log<-cbind(n_min_dist_log, mean_min_dist_log, sd_min_dist_log, min_min_dist_log, max_min_dist_log)
taba4<-rbind(all_rol_inverse, all_trade_total_log, all_alliance, all_refugees_log, all_min_dist_log, all_democracy, all_state_capac, all_gdppc_log, all_pop_log, all_right_orient, all_civ_lib, all_orign_country_past, all_press_freedom)
taba4<-as.data.frame(taba4)
names(taba4) <- gsub("_rol_inverse", "", names(taba4))
taba4<-round(taba4, digits=2)
names<-c("Rule of Law (inverse)", "Total Trade (log)", "Alliances", "Refugees (log)", "Minimum Distance (log)", "Democracy", "State Capacity", "GDPPC (log)", "Population (log)", "Right Orientation", "Civil Liberties", "Origin Country", "Press Freedom")
taba4<-cbind(names,taba4)

# Print Descriptive Stats 
taba4

#############################################################################################################################
# Table A.5: Zero-inflated negative binomial regression, interaction between the rule of law (inverse) and total trade (log)
#############################################################################################################################

# Model 1: Baseline model, independent variables with interaction
taba5_m1 <- zeroinfl(tr_count_det_ren_dep ~ rol_inverse*trade_total_log + rol_inverse + trade_total_log | 1, dist = 'negbin', data = tr_rd)

# Model 2: Full model, independent variables with interaction and control variables
taba5_m2 <- zeroinfl(tr_count_det_ren_dep ~ rol_inverse*trade_total_log + rol_inverse + trade_total_log + alliance + refugees_log + min_dist_log + democracy + state_capac + gdppc_log + pop_log + right_orient + civ_lib + orign_country_past + press_freedom | 1, dist = 'negbin', data = tr_rd)

# Print Regression Table
screenreg(list(taba5_m1, taba5_m2), custom.model.names = c("Model 1", "Model 2"), stars = c(0.01, 0.05, 0.1), digits = 3)

# Observations
taba5_n_m1<-taba5_m1$n
taba5_n_m2<-taba5_m2$n

# Log Likelihood
taba5_loglik_m1<-taba5_m1$loglik
taba5_loglik_m2<-taba5_m2$loglik

# theta 
taba5_theta_m1<-taba5_m1$theta
taba5_theta_m2<-taba5_m2$theta

# Prob. > X2
taba5_px2_m1<-lrtest(taba5_m1)
taba5_px2_m1<-taba5_px2_m1$`Pr(>Chisq)`[2]
taba5_px2_m2<-lrtest(taba5_m2)
taba5_px2_m2<-taba5_px2_m2$`Pr(>Chisq)`[2]

# Create Model Fit Stats Table
taba5_m1_fit<-rbind(taba5_n_m1, taba5_loglik_m1, taba5_theta_m1, taba5_px2_m1)
taba5_m2_fit<-rbind(taba5_n_m2, taba5_loglik_m2, taba5_theta_m2, taba5_px2_m2)
taba5_fit<-round(cbind(taba5_m1_fit,taba5_m2_fit), digits=3)

# Print Model Fit Stats
taba5_fit

################################################################################################################################################################################################################################
# Table A.6: Zero-inflated negative binomial regression, count of transnational repression incidences, excluding instances where it is less clear that the origin country government was ultimately responsible for the incident
################################################################################################################################################################################################################################

# Model 1: Baseline model, independent variables
taba6_m1 <- zeroinfl(tr_count_political ~ rol_inverse + trade_total_log | 1, dist = 'negbin', data = tr_rd)

# Model 2: Full model, independent variables and control variables
taba6_m2 <- zeroinfl(tr_count_political ~ rol_inverse + trade_total_log + alliance + refugees_log + min_dist_log + democracy + state_capac + gdppc_log + pop_log + right_orient + civ_lib + orign_country_past + press_freedom | 1, dist = 'negbin', data = tr_rd)

# Print Regression Table
screenreg(list(taba6_m1, taba6_m2), custom.model.names = c("Model 1", "Model 2"), stars = c(0.01, 0.05, 0.1), digits = 3)

# Observations
taba6_n_m1<-taba6_m1$n
taba6_n_m2<-taba6_m2$n

# Log Likelihood
taba6_loglik_m1<-taba6_m1$loglik
taba6_loglik_m2<-taba6_m2$loglik

# theta 
taba6_theta_m1<-taba6_m1$theta
taba6_theta_m2<-taba6_m2$theta

# Prob. > X2
taba6_px2_m1<-lrtest(taba6_m1)
taba6_px2_m1<-taba6_px2_m1$`Pr(>Chisq)`[2]
taba6_px2_m2<-lrtest(taba6_m2)
taba6_px2_m2<-taba6_px2_m2$`Pr(>Chisq)`[2]

# Create Model Fit Stats Table
taba6_m1_fit<-rbind(taba6_n_m1, taba6_loglik_m1, taba6_theta_m1, taba6_px2_m1)
taba6_m2_fit<-rbind(taba6_n_m2, taba6_loglik_m2, taba6_theta_m2, taba6_px2_m2)
taba6_fit<-round(cbind(taba6_m1_fit,taba6_m2_fit), digits=3)

# Print Model Fit Stats
taba6_fit

#######################################################################################################################################################################################
# Table A.7: Zero-inflated negative binomial regression, count of transnational repression incidences, excluding instances where individuals are accused of terrorism and/or extremism
#######################################################################################################################################################################################

# Model 1: Baseline model, independent variables
taba7_m1 <- zeroinfl(tr_count_terror ~ rol_inverse + trade_total_log | 1, dist = 'negbin', data = tr_rd)

# Model 2: Full model, independent variables and control variables
taba7_m2 <- zeroinfl(tr_count_terror ~ rol_inverse + trade_total_log + alliance + refugees_log + min_dist_log + democracy + state_capac + gdppc_log + pop_log + right_orient + civ_lib + orign_country_past + press_freedom | 1, dist = 'negbin', data = tr_rd)

# Print Regression Table
screenreg(list(taba7_m1, taba7_m2), custom.model.names = c("Model 1", "Model 2"), stars = c(0.01, 0.05, 0.1), digits = 3)

# Observations
taba7_n_m1<-taba7_m1$n
taba7_n_m2<-taba7_m2$n

# Log Likelihood
taba7_loglik_m1<-taba7_m1$loglik
taba7_loglik_m2<-taba7_m2$loglik

# theta 
taba7_theta_m1<-taba7_m1$theta
taba7_theta_m2<-taba7_m2$theta

# Prob. > X2
taba7_px2_m1<-lrtest(taba7_m1)
taba7_px2_m1<-taba7_px2_m1$`Pr(>Chisq)`[2]
taba7_px2_m2<-lrtest(taba7_m2)
taba7_px2_m2<-taba7_px2_m2$`Pr(>Chisq)`[2]

# Create Model Fit Stats Table
taba7_m1_fit<-rbind(taba7_n_m1, taba7_loglik_m1, taba7_theta_m1, taba7_px2_m1)
taba7_m2_fit<-rbind(taba7_n_m2, taba7_loglik_m2, taba7_theta_m2, taba7_px2_m2)
taba7_fit<-round(cbind(taba7_m1_fit,taba7_m2_fit), digits=3)

# Print Model Fit Stats
taba7_fit

########################################################################################################################################
# Table A.8: Zero-inflated negative binomial regression, count of transnational repression incidences, including Interpol abuse variable
########################################################################################################################################

# Model 1: Baseline model, independent variables
taba8_m1 <- zeroinfl(tr_count_det_ren_dep ~ tr_count_interpol + rol_inverse + trade_total_log | 1, dist = 'negbin', data = tr_rd)

# Model 2: Full model, independent variables and control variables
taba8_m2 <- zeroinfl(tr_count_det_ren_dep ~ tr_count_interpol + rol_inverse + trade_total_log + alliance + refugees_log + min_dist_log + democracy + state_capac + gdppc_log + pop_log + right_orient + civ_lib + orign_country_past + press_freedom | 1, dist = 'negbin', data = tr_rd)

# Print Regression Table
screenreg(list(taba8_m1, taba8_m2), custom.model.names = c("Model 1", "Model 2"), stars = c(0.01, 0.05, 0.1), digits = 3)

# Observations
taba8_n_m1<-taba8_m1$n
taba8_n_m2<-taba8_m2$n

# Log Likelihood
taba8_loglik_m1<-taba8_m1$loglik
taba8_loglik_m2<-taba8_m2$loglik

# theta 
taba8_theta_m1<-taba8_m1$theta
taba8_theta_m2<-taba8_m2$theta

# Prob. > X2
taba8_px2_m1<-lrtest(taba8_m1)
taba8_px2_m1<-taba8_px2_m1$`Pr(>Chisq)`[2]
taba8_px2_m2<-lrtest(taba8_m2)
taba8_px2_m2<-taba8_px2_m2$`Pr(>Chisq)`[2]

# Create Model Fit Stats Table
taba8_m1_fit<-rbind(taba8_n_m1, taba8_loglik_m1, taba8_theta_m1, taba8_px2_m1)
taba8_m2_fit<-rbind(taba8_n_m2, taba8_loglik_m2, taba8_theta_m2, taba8_px2_m2)
taba8_fit<-round(cbind(taba8_m1_fit,taba8_m2_fit), digits=3)

# Print Model Fit Stats
taba8_fit

###################################################################################################################################################################
# Table A.9: Zero-inflated negative binomial regression, count of transnational repression incidences, excluding Thailand (host country) and China (origin country)
###################################################################################################################################################################

# Reduce to sample without Thailand and China dyads
taba9_df<-tr_rd[tr_rd$ccode1!=800 & tr_rd$ccode2!=710,]

# Model 1
taba9_m1 <- zeroinfl(tr_count_det_ren_dep ~ rol_inverse + trade_total_log | 1, dist = 'negbin', data = taba9_df)

# Model 2
taba9_m2 <- zeroinfl(tr_count_det_ren_dep ~ rol_inverse + trade_total_log + alliance + refugees_log + min_dist_log + democracy + state_capac + gdppc_log + pop_log + right_orient + civ_lib + orign_country_past + press_freedom | 1, dist = 'negbin', data = taba9_df)

# Print Regression Table
screenreg(list(taba9_m1, taba9_m2), custom.model.names = c("Model 1", "Model 2"), stars = c(0.01, 0.05, 0.1), digits = 3)

# Observations
taba9_n_m1<-taba9_m1$n
taba9_n_m2<-taba9_m2$n

# Log Likelihood
taba9_loglik_m1<-taba9_m1$loglik
taba9_loglik_m2<-taba9_m2$loglik

# theta 
taba9_theta_m1<-taba9_m1$theta
taba9_theta_m2<-taba9_m2$theta

# Prob. > X2
taba9_px2_m1<-lrtest(taba9_m1)
taba9_px2_m1<-taba9_px2_m1$`Pr(>Chisq)`[2]
taba9_px2_m2<-lrtest(taba9_m2)
taba9_px2_m2<-taba9_px2_m2$`Pr(>Chisq)`[2]

# Create Model Fit Stats Table
taba9_m1_fit<-rbind(taba9_n_m1, taba9_loglik_m1, taba9_theta_m1, taba9_px2_m1)
taba9_m2_fit<-rbind(taba9_n_m2, taba9_loglik_m2, taba9_theta_m2, taba9_px2_m2)
taba9_fit<-round(cbind(taba9_m1_fit,taba9_m2_fit), digits=3)

# Print Model Fit Stats
taba9_fit